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Abstract 

A method to incorporate passive microwave remote sensing measurements within a 
spatially distributed snow hydrology model to provide estimates of the spatial distribution 
of snow water equivalent (SWE) as a function of time is implemented. The passive 
microwave remote sensing measurements are at 25 km resolution. However, in mountain 
regions the spatial variability of SWE over a 25 km footprint is large due to topographic 
influences. On the other hand, the snow hydrology model has built-in topographic 
information and the capability to estimate SWE at a 1 km resolution. In our work, the 
snow hydrology SWE estimates are updated and corrected using SSM/I passive microwave 
remote sensing measurements. The method is applied to the Upper Rio Grande River 
Basin in the mountains of Colorado. The change in prediction of SWE from hydrology 
modeling with and without updating is compared with measurements from two SNOTEL 
sites in and near the basin. The results indicate that the method incorporating the remote 
sensing measurements into the hydrology model is able to more closely estimate the 
temporal evolution of the measured values of SWE as a function of time. 


1. Introduction 

Mountain seasonal snowpacks provide much of the potable water supply in many 
areas, particularly the western United States. These mountain snowpacks exhibit high 
spatial variability in their properties, especially snow water equivalent (SWE), which is an 
important variable in hydrology and water resources engineering. The total water content 
of the snow and its spatial distribution determine the amount and timing of snow meltwater 
available for water supply and hydropower generation. The development of a method to 
estimate spatial distributions of snow properties such as SWE at a 1 km spatial resolution 
will provide information useful to snow melt prediction. 

The estimation of distributed SWE, however, remains a challenge due to the high 
variability of this quantity over small distances. Point measurements of snow properties, 
such as the SNOTEL (SNOwpack TELemetry) measurements of SWE are sparsely 
distributed and not representative of the parameter spatial variability. This makes it difficult 
to obtain a sufficiently accurate estimate of the water volume available in the snowpack and 
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to incorporate these measurements into basin-scale spatially distributed hydrology models 
which operate on a spatial grid of resolution 1 km or less. 

Satellite remote sensing appears to have the potential for overcoming this deficiency 
due to its inherent ability to measure spatial averages of snow properties. By their nature, 
satellite remote sensing measurements can provide spatial estimates of snow properties and 
are suited as inputs to spatially distributed hydrologic models. There has been previous 
success in using satellite remote sensing measurements for the determination of snow 
properties. The most successful has been the delineation of snow covered area (SCA) 
using either visible (Dozier and Marks, 1987) or microwave (Josberger and Beauvillain, 
1989) satellite measurements. The estimation of liquid water in the snowpack using 
microwave radar measurements has also shown some success when compared with field 
measurements (Shi and Dozier, 1995). The capabilities of satellite derived SCA within a 
snow hydrology model has been shown (i.e., Rango, 1988). More information on snow 
remote sensing and its applications within snow hydrology models can be found in the 
recent reviews of Rango (1993), Engman (1995), and Bales and Harrington (1995). 

Passive microwave remote sensing should also be able to provide information on 
other snow properties. The received microwave signal from the earth's surface is 
dependent on snow properties such as grain size, temperature, depth, and density. Unlike 
the use of SCA within hydrology models, less has been done in fully exploiting this 
information with a snow hydrology model. Rango et al. (1989) and Chang et al. (1991) 
have explored the use of passive microwave remote sensing measurements to obtain basin- 
wide estimates of SWE. Examinations of the relation between snow properties and passive 
microwave brightness temperatures have been undertaken (see for example Hallikainen and 
Jolma, 1992; Wang et al., 1992), but the relationship is not direct and estimation of snow 
properties remains difficult. One of the difficulties in implementing a direct estimation 
technique for snow properties is the many-to-one inverse problem encountered in going 
from the measured microwave brightness temperatures to the snowpack properties. 
Another problem is the difference in spatial resolution between a spatially distributed 
hydrology model (1 km, or less) and the footprint resolution of the passive microwave 
measurements (25 km). Thus, these previous algorithms retrieve averaged snow 
parameters over the 25 km footprint. 

To tackle these problems we use a method that combines several elements: a 
spatially distributed snow hydrology model, a snow-microwave radiative transfer model, 
and a parameter estimation technique with error minimizing steepest descent gradient 
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search. The spatial hydrology model operates at a 1 km pixel resolution using point 
meteorological measurements as forcings and outputs snowpack parameters at this 
resolution. Based on the snow property outputs, a microwave scattering model is used to 
determine the brightness temperatures for each pixel. The 1 km pixel estimates are 
integrated to the footprint measurement resolution and compared with the measured SSM/I 
brightness temperatures. The adjusted brightness temperatures are disaggregated to the 
pixel level and inverted to give new snow parameters used as updating inputs to the snow 
hydrology model. The hydrology model continues running until the next update. 

To test this method, we apply it to a hydrological basin in the mountains of 
Colorado. The estimates of the temporal evolution of SWE are compared to SNOTEL 
ground truth measurements at two sites. 

2. Methodology 

The methodology outlined in the introduction section consists of combining several 
sources of information to provide the best estimate of the snowpack state parameters as 
possible. We use a spatially distributed snow hydrology model and microwave satellite 
measurements. These are linked by a snow-microwave scattering model and a parameter 
estimation technique. In this section the new algorithm is outlined with each component 
discussed in greater detail in the next section. The procedure used in obtaining the spatial 
snow water equivalent estimates is given in the flow chart of Figure 1. Each step is 
denoted by a letter in the flow chart and corresponds to the lettered description below. 

a) The spatial data files are input to the snow hydrology model to determine the pixel 
characteristics of the model grid at 1 km resolution. The hourly point meteorological 
station information is used to drive the energy balance component of the snow hydrology 
model. These measurements consist of air temperature, dew point temperature, surface 
pressure, precipitation, cloud cover, and wind speed. These forcing variables are needed 
as input for each pixel and are adjusted for each pixel based on its topographic 
characteristics. 

b) Net energy fluxes are calculated for each pixel using the snowpack surface energy 
balance model which operates at a 1 km spatial resolution. The output from the snow 
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hydrology model is snow depth, density, grain size, and temperature. These parameters 
are determined for each pixel on an hourly basis. 

c) At the end of each day, the four output snowpack parameters are saved for use as input 
to the remote sensing-electromagnetic scattering model. The end-of-day values coincide 
with the nighttime satellite measurements used. 

d) Based on the snowpack parameters obtained from the snow hydrology model, the 
brightness temperature estimates at 1 km spatial resolution are made. This is done for 19 
and 37 GHz for both vertical and horizontal polarizations. Each is integrated to the 25 km 
scale to match the scale of the satellite measurements. 

e) SSM/I measurements, with a resolution of 25 km, are used to compare with those 
obtained from the models. The area of coverage is determined for each remote sensing 
measurement; this is the overlap between the spatial domain of the models and the field of 
view of the sensor measurements. 

0 The differences between the spatially averaged simulated brightness temperatures and 
those measured by the sensor are determined for each coverage area. The differences are 
uniformly distributed back to each individual pixel in the respective coverage area. Each 
pixel is adjusted by the same difference so that the average value of the pixel brightness 
temperatures for the coverage area equals the values obtained from the SSM/I 
measurements. 

g) Snow parameters at the pixel (1 km) scale are determined from the adjusted brightness 
temperatures through inversion. The updated snow parameters are returned to the snow 
hydrology model which moves forward in time until the next available nighttime sensor 
observation, whereupon the process repeats itself. 

Our algorithm is a systematic combination of hydrology model and remote sensing 
measurements. The remote sensing measurements are used to update and correct the 
hydrology model prediction. On the other hand, the hydrology model constrains the 
parameter inversion from the remote sensing measurements. This solves the many-to-one 
inverse problem in remote sensing, as the hydrology model gives a priori estimates of the 
snow parameter inversion. In addition, this algorithm gives snow parameter products at 1 
km resolution as well as the temporal evolution of snow parameters. 

The available data is discussed in the next section. Due to a lack of ground truth, 
we implement this algorithm only for selected points where measurements of snow water 
equivalent are available. This is used to validate the algorithm and demonstrate that the 
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algorithm has the potential to provide snow parameter estiates at the 1 km pixel resolution. 
Thus, we have not implemented the 25 km integration of step (e) and all of step (f). 


3. Data Characteristics 

The basin chosen for the application of the technique is the upper portion of the Rio 
Grande River in the mountains of Colorado. The basin is defined by the streamgauge at 
Del Norte, Colorado and has an area of 3419 km 2 with an elevation range of 2432-4215 
m. This basin has been used for previous snow-remote sensing studies (e.g., Rango et ah, 
1989; Chang et ah, 1991). The location of the basin is shown in Figure 2. We use the 
1992-1993 snow year for the application. We employ several data types in the 
implementation and validation of the new algorithm: hydrology inputs, remote sensing 
measurements, and ground truth. These are discussed below and listed in Table 1. 

A. Hydrology model inputs 

The data used to run the hydrology model consists of topographic and weather data. 
The basin topographical characteristics are determined from 30 arc-second digital elevation 
model (DEM) output, which gives an approximate pixel resolution of 1 km and defines the 
spatial grid of the hydrology model. The basin area and its topographic representation are 
determined from the DEM and shown in Figure 3. From the DEM file, pixel elevation is 
determined, and slope, aspect, and shading characteristics of the pixels are calculated. 

The forest cover plays an important role in the energy and mass balances of the 
snowpack through the reduction of shortwave radiation and wind speed, and the 
interception of snowfall by the forest canopy. Thus it is important to include forest cover 
information in the snow hydrology model. Pixel forest cover was determined from the 
Land Use and Land Cover digital files from the U.S. Geological Survey (USGS, 1990), 
which uses a classification scheme developed by Anderson et al. (1976). The land use 
pixels, with 200 m resolution, are aggregated to the DEM pixel scale (1 km) and used to 
determine whether each hydrology model pixel is open, forested, or some fraction thereof. 

The snow hydrology model uses meteorological inputs to determine the 
precipitation characteristics as well as drive the snowpack surface energy balance. Hourly 
point meteorological measurements from the National Weather Service are used in 
conjunction with a distribution scheme based on the topography to provide spatial inputs to 
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the snow hydrology model. These consist of air temperature, wind speed, relative 
humidity, cloud cover, and precipitation. Information from additional sites that record 
daily precipitation totals is also used. Additionally, U.S. Geological Survey daily 
streamflow values at the outlet of the basin are used to estimate the volume of precipitation 
over the basin. 

B. Satellite remote sensing measurements 

Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave 
Imager (SSM/I) brightness temperatures are obtained from the Distributed Active Archive 
Center (DAAC) at the Marshall Space Flight Center. The SSM/I measures microwave 
emission at four frequencies: 19, 22, 37 and 85 GHz with dual-polarization, except at 22 
GHz, which is vertical only. More information on this sensor can be found in Hollinger et 
al. (1990). We selected the F8 satellite, using the brightness temperature measurements at 
the 19 and 37 GHz frequencies with horizontal and vertical polarizations, providing four 
different measurements for each footprint. Each day with measurements over the basin is 
determined. The nighttime measurements are saved and the daytime measurements are 
discarded due to the possibility of small amounts of liquid water present in the snowpack as 
a result of slight surface melting during the day, which re-freezes at night. 

While the orbit of the satellite is sun-synchronous, providing approximately the 
same time of coverage over the basin, there are some days for which there is no satellite 
coverage of the basin. Additionally the footprint positions vary slightly from day to day. 
This means that the relationship between the satellite footprint that covers a pixel and the 
pixel snow properties is not constant since the distribution of pixels within the footprint 
changes. The change in footprint location for a SNOTEL pixel (#4) over four consecutive 
days is shown in Figure 4. This effect must be taken into account when attempting to relate 
pixel values of snow properties with the brightness temperature measurements of the 
satellite footprint and will be discussed further in a later section. 

C. Ground truth measurements for algorithm validation 

The SNOTEL measurement network records daily SWE values at selected points 
throughout the mountains of the western United States. Six of these stations are located 
within or very near the Upper Rio Grande basin for the year of study. Their names, 
locations, and elevations are listed in Table 2 and their locations with respect to the basin 
shown in Figure 3 (as indicated by an 'X'). The stations are all in the upper area of the 
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basin where the snow amounts are larger and measurements more important. For this 
study, the SNOTEL values are taken as representing the SWE for the 1 km pixel within 
which it lies. These measurements represent the ground truth to evaluate the performance 
of and validate the algorithm. 


4. Implementation of Methodology 

The basic idea of combining snow hydrology information with passive microwave 
remote sensing measurements through parameter estimation techniques was outlined in a 
previous section. There are several ways of implementing the three major components of 
this algorithm. Several different types of snow hydrology models and microwave radiative 
transfer models are available. Additionally there are several possible approaches to the 
parameter estimation procedure of error minimization. We have chosen the following ways 
to implement the three components previously discussed. 

A. Snow hydrology model 

The snow model contains several components: accumulation, snow surface energy 
balance, and internal snowpack physics. Snowmelt runoff is not currently considered 
since we investigate only the accumulation portion of the snow season. The model 
operates on a spatial grid of resolution approximately 1 km as determined by the DEM. 
Hourly meteorological inputs are input to the snow hydrology model to drive the 
accumulation and surface energy balance components. These inputs are adjusted for each 
model pixel using an elevation relationship such as air temperature lapse rate, 
dT a /dz = constant. 

The variables describing the snowpack for use in the radiative transfer model are 
snow depth d , snowpack temperature 7j, grain size (radius) g x , and fractional volume f s . 
These four variables are referred to as the snowpack state variables S. Each pixel within 
the spatially distributed snow hydrology model will have its own set S. The snow 
hydrology model maintains two related snowpack state variables: snow density p s and 
snow water equivalent IV. These are related to the parameters in S by 

ft =0.91/,, (1) 

and 
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W = p s d, (2) 

where 0.9 1 is the assumed density of ice ( g • cm' 5 ) and W is in cm. These state variables 
evolve over time as a function of the mass and energy balances as discussed below. 

The model operates over a discrete time step At, during which the snowpack water 
equivalent evolves according to 

W(t + At) = W(t)+P-M-E, (3) 

where P is snowfall, M is melt, and E is the loss due to sublimation, all in units of cm. 

The accumulation component of the model, P, is based on point precipitation and 
temperature records in the area, to determine the amount that falls and whether it falls as 
rain or snow. The precipitation/elevation relationship is adjusted using the SNOTEL data 
by adjusting the total volume of precipitation at the point. The total annual volume of 
precipitation over the basin is adjusted using the streamflow measurements. 

Following Marks (1988) we use a two-layer snowpack model, which consists of an 
"active layer" from which melt is produced, and a secondary layer which contains the 
remainder of the snow. The active layer is chosen to have a depth of 10 cm, thus capturing 
most of the depth of solar radiation penetration (approximately 90 percent), while not 
requiring an unduly large amount of heat to initiate melt conditions. The melt from the 
active layer is 


M = — — , (4) 

An 

where L m = 79.7 cal-g -1 is the latent heat of melt (the amount of energy required to 
transform ice to water), Q n is the net energy flux into the snowpack, and Q cc is the cold 
content of the active layer, representing the amount of heat necessary to raise its 
temperature to 0°C. If Q cc > Q n , then no melt occurs. Since the model only considers the 
accumulation and ablation process, the meltwater is immediately removed from the 
snowpack without any routing, and the possibility of the meltwater refreezing in the pack is 
not considered. 

The snowpack cold content is 
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Qcc = ~C S WT S , (5) 

where C s = 0.5 cal g' 1 deg' 1 is the snowpack specific heat. If the net energy into the 
snowpack is great enough to raise the snowpack temperature to 0°C {Q n >Q cc ) then melt 
will be produced according to (4). If no melt is produced, the updated snowpack 
temperature is 


T, (• + &) = + (6) 

The mass loss due to sublimation is calculated as 



where Q e is the latent heat transfer from the energy balance and L ( = 677 cal g' 1 is the 
latent heat of sublimation (the amount of energy required to transform ice to vapor). 

The net energy flux at the surface, Q n , is given by 

Q„ =Qt + Q< + Q h + Q t , ( 8 ) 

where the fluxes Q s , Q n Q h , Q r are respectively shortwave, longwave, sensible, and latent 
heat. The ground-to-snowpack and the precipitation-to-snowpack heat transfers are small 
relative to the other terms of (8) and are thus omitted. 

The equations used to estimate the components of the snow surface energy balance 
are based on the point model from Marks (1988). However, it has been generalized to a 
spatial model incorporating the topography and forest cover characteristics of the basin. 
This includes accounting for the topographical and vegetative effects on solar and longwave 
radiation. More information on the computation of the different components in the surface 
energy balance (Equation 8) is in Appendix A. SWE outputs are compared with the 
SNOTEL stations to determine that they are reasonably close. 

The snowpack physics algorithms, which describe the evolution of snowpack 
properties such as snow grain size and density are taken from the SNTHERM model of 
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Jordan (1991). The snow metamorphism is driven by the outputs from the energy and 
mass balances previously discussed. The algorithms are simplified to operate for a one- 
layer snowpack; there is no variability in the vertical direction. 

The equation for the evolution of snow grain size is 

g,(r + A/) = 0.14-^— Ar + g,(r), (9) 

4 *,(0 

where g 2 is an empirical constant with a value of 4.0x l(T 8 m 2 s' 1 . 

There are two components in the snow compaction process: destructive 

metamorphism and overburden. Destructive metamorphism is important initially when 
snow density is less than 0.15 g em' 3 . For dry snow, the empirical function used for 
compaction is 


1 dAz 


A z dt 


m 


= -2.778 x 10^ xq xc 4 x exp[-0. 04(273. 15 -T)], 


( 10 ) 


where c 4 is 1 for dry snow (the case here) and 


c 3 = 1 

c 3 =exp[-46.(p i -0.15)] 


if p s <0.15 g em' 3 , 

if p s >0.15 g em' 3 . 


After snow has undergone its initial settling stage, densification proceeds at a slower rate 
determined mainly by compaction due to snow overburden pressure P s such that 


1 dAz 
A z dt 0 


- exp[-c 5 (273. 1 5 - 7)]exp(-c 6 • p , ) , 
»0 


( 11 ) 


where rj 0 is a viscosity coefficient for snow and c 5 and c 6 are constants based on 
observational evidence. The compaction rate is thus the combination of these two terms. 
As the depth changes but the water equivalent does not (no melt), the density of the 
snowpack will also change. 

The final output of four snowpack parameters for the microwave scattering model 
is: d snow depth, T s snowpack average temperature, g t snowpack average grain size, and 
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/., snowpack average fractional volume, S = {d, T s ,g t ,f s }. These four snowpack 
parameters serve as the link between the microwave scattering model and the snow 
hydrology model. 

Figure 5 shows the comparison between the hydrology model output SWE and the 
measured SWE for two SNOTEL stations: #2 which is a high elevation station and #4 
which is a low elevation station. The model is able to reproduce the snowfall events, but is 
usually not able to accurately reproduce the amounts. This is due to the simple scheme 
used to determine snowfall amounts in a topographically complex area. This is where the 
incorporation of information from passive microwave remote sensing measurements can 
improve SWE estimates. 

To demonstrate the evolution of these parameters in time, simulation outputs from 
the hydrology model for all four snow parameters at station #2 (Figure 6a-d) and station #4 
(Figure 7a-d) are shown. In both cases the depth increases whenever there is a snowfall 
event, then decreases as the snow compacts. Snowpack temperature stays low until the 
end of March, when it starts to increase as the melt season approaches. Grain radius and 
fractional volume both follow nearly the same shape of curve. This is due to the averaging 
of the snowpack properties to one layer. These values increase as the snowpack 
metamorphoses, until a new snowfall event, which reduces the average snowpack value by 
introducing fresh snow, which has smaller grain sizes and is less dense. 

B. Snow-microwave radiative transfer 

In passive microwave remote sensing the microwave emission from the ground 
surface is measured. This radiation is attenuated by the snowpack as a function of its 
temperature, depth, density, and grain size, as well as liquid water. To avoid the 
complications due to liquid water in the snowpack we use only nighttime measurements, 
when the snowpack has re-frozen, and do not continue the snow simulation into the melt 
season. From the snow hydrology model of the previous section, we have estimates of the 
four snowpack parameters S at the pixel scale. The brightness temperatures at this 1 km 
scale can be estimated using a microwave scattering model. Thus, from the model we have 
the brightness temperature vector T b = \T h ^ , T bnit , T hm , T h> ^ J, where 19 and 37 
represent the frequency in GHz and V and H the vertical and horizontal polarization. 

The model used to simulate the microwave interaction with the snowpack is the 
dense medium radiative transfer (DMRT) model (Tsang 1987, 1992). Unlike the 
traditional radiative transfer models, the DMRT takes into account the dependency of 
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scattering upon relative particle positions. This is important in a dense medium such as 
snow. The model incorporates a particle size distribution using a modified gamma 
distribution. Thus a medium containing particles of different sizes, such as snow, can be 
accounted for. The model accounts for scattering by particles as well as the interaction at 
the snow/soil interface and the snow/air interface. The equations describing the model are 
presented below. Additional details on the model and its implementation are in Tsang et al. 
(1992). 

We consider thermal emission from a layer of dielectric particles of permittivity e s 
embedded in a background medium of permittivity e overlaying a homogeneous half-space 
of dielectric medium of permittivity e 2 (see Figure 8). Note for the case of snow, the 
background permittivity e is e 0 of free space and the particle permittivity e t will be that of 
ice. The particle sizes obey a size distribution n(g) which is the number of particles per 
unit volume with radii between g and g + dg. The medium is of uniform temperature T. 
Then the dense media radiative equation for passive remote sensing assumes the following 
matrix form, for 0 < 9 < k\ 

d - - rr 

cos 9 — I{z,9) = -kJ(z,6)+ k [\-(b)CT 
oz I_1 

K'd) Jjsin 9V(9 , 9') • /(z, 9')dd' (12) 


where 


Hz,9) 


K(z,9) 

I„(z,d) 


(13) 


and I v , I h are the vertical and horizontal specific intensities, respectively. Also in (12), 
C=K,K' , /(XV), K h is Boltzman's constant, A is the free space wavelength, K' is the 
real part of the effective wave number in region 1, it is the free space wave number, 0) is 
the albedo, and K e = 2 Im( K) is the extinction rate of the specific intensity. In ( 12) 


P (9,9') 


Pu{e,9') P i2 {9,9') 

Pv(M') P2 *{W) 


(14) 


where 


p u (9,9') = 2 sin 2 0sin 2 9' + cos 2 Ocos 2 9' 


(15) 
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Pi 2 {0,9') = cos 2 9 (16) 

p 2l (d,6') = cos 2 d' ( 17 ) 

p 22 (0,0') = l. (18) 

For 0 < 6 < k\ 2 , the boundary conditions for (12) at the air-snow interface, z = 0, are 

I v (z = 0,K-d) = R v (d)I v {z = 0,e) (19) 

and 

i h (z = o,7r-e) = R h (e)i h (z=o,9), ( 20 ) 

and for the snow-ground interface, z = -d, they are 

iXz = -d,e) = R Vft {e)i v {z = -d,K-d)+(i-^{d))cT g (2i) 

and 

h (z = ~d, 6) = R hg (e)I h (z = -d,x-6) + ( I - R hx (0))CT g , ( 22 ) 

where R v , R h , ^ e ,and R hg are reflectivities. 


After (12) is solved subject to the boundary conditions of ( 1 9)— (22), the brightness 
temperatures in the direction d Q , where B 0 = sin'^AT'sin d/k) is related to 0 by Snell’s 
law, in region 0 for vertical and horizontal polarizations are given by 


r 

o 

CO 

> 

CQ 

1 

_ i 

'( 1 - *,(»))/,(* - 0 . 9 )' 

1 

o 

i 

~ c 

_(l-R ft (9)K(z = O,0)_ 


(23) 


The differences between the dense medium theory and the conventional radiative 
transfer theory are the calculations of K, the extinction rate K e , and the albedo To in terms 
of the physical parameters of the medium, which are represented by e s and the size 
distribution n(g). These parameters are determined by using approximations of the Dyson 
and Bethe-Salpeter equations. The cross pair distribution functions of multiple particle 
sizes are calculated through the Percus-Yevick approximation (Percus and Yevick, 1958; 



15 


Baxter, 1970) that expresses the correlations of particle positions in terms of the size 
distribution n(g). 

Also it is necessary to consider the effects of surface roughness in the Fresnel 
reflectivities, especially when the snowpack is thin. Following Wang et al. (1983) we use 
two polarization mixing parameters, one for each interface. For the upper boundary, the 
air-snow interface, z = 0, 


R.(0)=[(i-ek(0)+Qr,(e)] 

(24) 

«,(e)=[(i-ek(e)+2r.(0)]. 

(25) 


where Q is the mixing fraction between the two reflection polarizations for the snow-air 
interface and r v and r h are the Fresnel reflectivities with effective propagation constant K 
determined by 


and 


r v (e) = 


r>(0) = 


k 2 cos 6 - K'[k 2 - K' 2 sin 2 df 


k 2 cos d + K'(k 2 - K' 2 sin 2 Of 
K'cos 6 - (k 2 - K' 2 sin 2 df 


K ' cos 6 + - K’ 2 sin 2 0) 


In the lower boundary, at the snow-ground interface, z = -d , 


(26) 


(27) 


and 


(28) 


with 


R ..( 0 )=[('-e,h( e ) + e,\( e )] 


(29) 
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and 



— k Q 2 cos 6 - K' 
e o 

( £ > 

-*-k 0 2 — K' 2 sin 2 9 

V e o ) 

X 

— k 0 2 cosd + K' 
£ o 

^k 2 -K' 2 sm-Q\ 

Uo J 

Ti 



K’ cose- 

( £ > 

— k 0 2 - K' 2 sin 2 6 

\ £ 0 ) 

¥> 

K' cosO + 

( £ 1 

^-k 2 -K' 2 sin 2 0 

V. £ o j 

IT ’ 


( 30 ) 


(31) 


where k Q is the free space wave number. 

In this study we use a Rayleigh size distribution of snow particles (West et al., 
1993) represented by 


n {g) = 


16( 5 ) 5 


exp 




(32) 


where ( g ) is the mean radius and / is the fractional volume of all the particles. The 
fractional volume and the mean radius are defined by 


f = C-f S 2 n{ S )d g 


(33) 


\;n{ g )dg- 


(34) 


The advantage of the Rayleigh size distribution is that there are only two parameters: / 
fractional volume, and (g) mean grain size. 

The medium input physical parameters of the dense medium radiative transfer 
equations using the Rayleigh size distribution include the previous vector of four snowpack 
parameters from the hydrology model S and a set of fixed parameters describing the 
electromagnetic properties of the dense medium, F, necessary to implement the DMRT. 
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This fixed set F consists of: 

1) the complex permittivity of the particles: e in = (3.2 + /0.00O25)e o and 

£ Ij7 = (3.2 + i0.001)e 0 , for 19 and 37 GHz respectively, 

2) the permittivity of the ground surface, e g , 

3) the ground surface temperature T g , 

4) the two previously described polarization reflectivity mixing parameters, Q and 

£2 k , and 

5) the maximum grain size in the grain size distribution of (32), g^. 

It is necessary to estimate some of the physical parameters of the scattering model, 
such as the emissivity of the soil and its initial temperature. This is done by calibration, 
using the microwave brightness temperature measurements for October, right before any 
snow has fallen and is covering the ground. This provides a calibration for soil properties, 
which affect the microwave brightness temperatures throughout the year, although much 
less as the snow season progresses and the snowpack becomes deeper. This allows for the 
estimation of £ g , T g , and the roughness mixing parameters, Q and Q g . 

Additionally, the determination of snow microwave properties through calibration is 
necessary. The maximum and minimum limits need to be set for the snow particle size 
distribution. The minimum value is based on the initial grain size used for fresh snow. 
The maximum value evolves through time, just as the shape of the particle distribution 
evolves through time. The hydrology model grain size (g) is taken as the distribution 
mean. The upper limit must be periodically re-established since as the mean grain size from 
the hydrology model evolves (Equation 9), the snow grain size distribution evolves with 
the upper limit increasing as well. The need to update is based on visual inspection of the 
snow grain size output from the model, see Figures 6c and 7c. The fixed parameter vector 
Fis {e 5 , e t , T g , Q , Q g , £_}. 

From these parameters, the brightness temperatures of vertical and horizontal 
polarization at specified frequencies and observation angles can be calculated from the 
above equations. Using the snow parameter output, the brightness temperatures at 19 and 
37 GHz for vertical and horizontal polarizations are obtained at the 1 km pixel scale. The 
efficient, approximate representation of the snow microwave radiative transfer model is 
described in Appendix B using a neural network representation. 
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C. Remote sensing data that overlap a pixel 

At the end of each day for which there is basin coverage from the sensor, the 
DMRT is used to estimate the four microwave brightness temperatures from the four 
snowpack parameters at pixel resolution (1 km). The pixel brightness temperatures are 
then averaged to the footprint scale (25 km) for each measurement. The difference between 
the pixel integrated brightness temperatures and SSM/I measured brightness temperatures is 
distributed back to the pixels from which these adjusted 1 km brightness temperatures can 
be used to obtain adjusted snowpack parameters. 

The critical point in the previously described procedure is in going from the 1 km 
scale to the 25 km scale and then back to the 1 km scale. The remote sensing 
measurements cover approximately a 25 km area and thus represent a spatial average of 
snow properties that may be quite heterogeneous. Depending on the location of the pixel 
with respect to the center of the footprint, the pixel-to-footprint measurement relationship 
can vary from day to day. 

To examine the possible fluctuations in remote sensing measurements due to 
footprint position, we compared SNOTEL sites with brightness temperatures assigned to a 
pixel based on the degree of coverage. The was examined as well in Wilson et al. (1996). 
Figure 9a shows the brightness temperatures for a coverage area with radius of 12.5 km, 
the actual footprint size. Figure 9b shows the brightness temperatures using a coverage 
radius of 4.0 km from the footprint center. As the pixel of interest becomes nearer to the 
center of the footprint, it should be more closely related to the footprint measurement. 
From the figures we see that the amount of fluctuation is reduced as the distance from pixel 
to footprint center is reduced, showing more clearly the expected trend of decreasing 
brightness temperatures with increasing snow depth. This method is employed in 
determining what measurements to use in the updating procedure. 


5. Results 

The implementation shown in Figure 1 is now demonstrated. A comparison 
between the two methods of determining SWE is made. The first method is the distributed 
snow hydrology model run independently of the remote sensing information, relying only 
on the topographic information and point meteorological measurements. The second 
method is using the snow hydrology model, but replacing the parameter outputs with the 
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updated set based on the remote sensing brightness temperatures. To determine how well 
the combined model performs compared to using only the spatially distributed snow 
hydrology model, we examine the results of the SWE curves from the two methods at two 
SNOTEL sites (#2 and #4). The comparison is made for the snow season of 1992-1993. 
The simulation starts on October 1, 1992, which is the first day of the hydrologic water 
year. The SNOTEL measurements are taken as the ground-truth for each site, representing 
the 1 km pixel within which it lies. 

The snow hydrology model is run first, independently of any remote sensing 
information, solely on the point meteorological forcings. Then the simulation is started 
with the remote sensing updating scheme. We start the updating when the SWE is greater 
than 5 cm so that the scattering of the microwave emission by the snow is significant. An 
update is made every five days. Linear interpolation is used between the brightness 
temperature observations that are used based on the previously discussed coverage 
criterion. In the time interval between updates, the two methods will follow each other 
since the SWE curves are dependent only on the precipitation estimates during the 
accumulation period. The methods will diverge only at the updating steps. 

The first station examined is station #4. It is a lower elevation station at 2860 m 
with smaller amounts of snow and maximum SWE of 15 cm. Its snow accumulation 
season begins later in the season than much of the basin. The results of the comparison are 
shown in Figure 10. We also include in the figure the measured SNOTEL SWE. One of 
the brightness temperature measurements, 37H, is included for reference. However, we 
have used all four brightness temperature measurements for the inversion of the snow 
parameters. Table 2 show the values for the fixed parameter vector F for this station. 

Since significant snow accumulation does not develop at this station until late 
December, the updating with remote sensing measurements is not started until December 
28. Improvement is obtained from using the updating procedure. Not every point is closer 
to the SNOTEL measurements than that obtained by the snow hydrology model alone, but 
more often it is closer. The updating is able to help compensate for poor snowfall estimates 
in the snow hydrology model, such as Day 1 10. The results in the figure indicate that the 
hydrology model with remote sensing updating is generally closer to the SNOTEL ground 
truth. 

Another test simulation is done for a high elevation station, #2, at 3420 m with 
larger amounts of snow and a SWE maximum of 40 cm. Note that a change in elevation of 
500 m between stations 2 and 4 gives a corresponding change in SWE of over 100%, 
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highlighting the spatial variability of SWE in this region. The comparison results are 
shown in Figure 11. Again an improvement is obtained from using the updating scheme. 
The results in the figure indicate that the hydrology model with remote sensing updating is 
generally closer to the SNOTEL ground truth. 

Between Day 70 and 75 a drastic drop is observed in measured brightness 
temperature. This is most probably due to the development of large grain sizes from a 
small melt event; the SNOTEL measurement drops temporarily during this time. This is 
one instance where it is important to recalibrate the DMRT, changing the upper grain size 
limit in the distribution. In general, every 15 days (after 3 updating events) the upper limit 
of the grain size g mM is recalibrated. Table 2 also shows values of the fixed parameter 
vector F for this site. The main difference is in the evolution of the grain size as evidenced 
by the differences in upper limit to the grain size distribution. The incorporation of the 
upper limit in the grain size distribution as a parameter in the neural network to avoid this 
recalibration is under research. 


6. Conclusions 

A method to combine passive microwave remote sensing measurements within a 
spatially distributed snow hydrology model to estimate the time evolution and spatial 
distribution of snow parameters has been presented. The method combines the information 
contained within the satellite microwave measurements with that obtained from the spatially 
distributed snow hydrology model. A spatially distributed snow hydrology model that 
uses spatial topographic inputs and point meteorological measurements to force the energy 
balance is combined with the SNTHERM model to obtain the four snowpack parameters 
for passive microwave remote sensing: snowpack density, depth, temperature, and grain 
size at a 1 km scale. The brightness temperatures for these 1 km snow parameters are 
estimated using the dense medium radiative transfer model. These are compared with the 
25 km resolution satellite measurements and the differences distributed back to the 1 km 
scale brightness temperatures. The adjusted snowpack parameters are determined by 
inverting the DMRT. 

The approach has enabled us to handle the two difficulties previously discussed: 
differences in spatial scale between the spatially distributed snow hydrology model (1 km) 
and remote sensing measurements (25 km), and the non-uniqueness of the inversion of the 
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remote sensing measurements. By implementing the 1 km resolution hydrology model we 
are able to relate the 1 km to the 25 km resolution SSM/I measurements with the DMRT 
and some estimation techniques. Inversion of the adjusted pixel brightness temperatures 
provides the adjusted snow parameters at the 1 km scale to update the snow hydrology state 
variables. Analysis of the satellite footprints is undertaken to determine effective pixel 
SWE relationships with degree of footprint coverage to that pixel. 

The approach is demonstrated for two SNOTEL sites, making a comparison 
between the case with only the snow hydrology model and the case that includes remote 
sensing updating. In each case there was overall improvement in the ability of the model to 
reproduce the SNOTEL SWE measurements when the remote sensing measurements were 
included. 

This methodology shows promising results in improving snow water equivalent 
amounts from snow hydrology models. This method can be easily extended to all pixels 
within the basin to map snow properties such as SWE. Additional remotely sensed types 
of information such as snow covered area can be incorporated within this framework. 
Bayesian updating techniques, which can provide better estimation are currently under 
investigation (Davis et al., 1995). 


Appendix A 

In this appendix the equations used to estimate the different components of the 
snow surface energy fluxes (Equation 8) are presented. These are for shortwave and 
longwave radiation, and latent and sensible heat. 

A. 1 Shortwave radiation 

Three components are considered in the calculation of shortwave radiation: direct, 
D, diffuse, 7, and sky reflected, R. Reflected radiation from neighboring slopes is 
neglected; an assumption that is justified for forested catchments, but would need to be 
retained in high alpine environments. The total shortwave radiation into the snowpack is 


Q$ = l + D + R 


(A. 1) 
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Direct radiation is that part of the solar beam that is not depleted due to scattering or 
absorption. Diffuse radiation is the fraction of the beam that is scattered in the forward 
direction. The sky reflected component is radiation reflected from the surface to the sky 
and then back to the surface again. This can be a significant component of total downward 
solar radiation over snow-covered surfaces where the albedo is high and there is low cloud 
cover. Presuming that measurements of /, D, and R are not available, which is normally 
the case, the direct and diffuse terms can be estimated using the clear sky equations of 
Munro and Young (1982) with the cloud effect estimated based on Davies et al. (1975) and 
Davies and Idso (1979). The direct component of solar radiation received at a snow 
covered pixel is 


/ = l 0 H { 1 - <*s )[ ^0 3 Vrs - (I 1 - Vw )] Vda Vds Vc F s (A.2) 

where I a is the direct solar radiation reaching the top of the earth’s atmosphere, H is a 
topographic parameter accounting for pixel slope, aspect, and shading, a s is the snow 
albedo, y is the transmissivity function for ozone absorption, Rayleigh scattering, water 
vapor absorption, aerosol absorption and scattering, and clouds (for subscripts O 3 , rs, w, 
da, ds, and C), and F s is an interception factor for forest cover. 

The diffuse radiation component is estimated assuming that the fraction of Rayleigh 
scattering in the forward direction is one half, and that the fraction of Mie scattering in the 
forward direction is 0.85 (Robinson, 1970). The diffuse radiation into the snowpack is 

0 ='o(l -“»){[ ^03 (1- «'r,)/2] + [v'Oj( 1 - V'r, ) - (1 - V,)] 

Vc V f F, (A3) 

where Vj- is the pixel view factor, which represents the fraction of the pixel open to the sky 
hemisphere. 

The reflected radiation is an infinite sum of reflected radiation between sky and 
atmosphere. However, it is only of significance in this study when there is cloud cover 
and only one return is considered. This is estimated as 


R = (I + D)a s {l-a s )a c V f 2 F s 2 


(A.4) 
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A. 2 Longwave radiation 

The longwave radiation term consists of four components: emissions from the 
atmosphere, clouds, forest, and snow surface. As with the shortwave component, the 
effect of the surface radiation from neighboring pixels is neglected, which is an appropriate 
assumption in the presence of forest cover. Longwave radiation is partitioned between 
forest-snow exchange and atmosphere-snow exchange. The longwave radiation balance at 
the surface is 


QL = {QLa+QLc)Vfi\-n + Q L fF-Qu (A. 5) 

where F is the fractional forest cover. The atmospheric component is 

Q[m ~ ^a < ^'a (A. 6a) 

the cloud component (Kimball et al., 1982) is 

n 

Qlc = 1L T s/8 e c, cr7 'c, 4 (A. 6b) 

(=1 

the forest component is 


Q Lf = £ f GT f 4 (A.6c) 

and the surface component is 

Qls = e s aT* (A.6d) 

where £ is emissivity, a is the Stefan-Boltzmann constant, T is temperature, and r 8 and 
/ 8 are the transmissivity and fraction of radiation in the 8-14 pm band, the region where 
the atmosphere is opaque in the infrared. The subscripts a , c,\/, and s represent air, cloud 
at level i, forest, and snow, respectively. 



24 


A.3 Heat transfer 

A bulk aerodynamic formulation (Sellers, 1965; Deardorff, 1968) is used for 
sensible heat: 


Qh - Pa c p^H u a^u{^a ~ T s )Sf (A.7) 

and for latent heat: 

fc = 0 ' 62 l P « L < C wU ,F„(e a - e s )S f (A.8) 

where C H and C w are the bulk transfer coefficients for heat and water, which depend on 
wind reference height, u a is the near-surface wind speed (e.g., at elevation 10 m above the 
surface), P a is the atmospheric pressure, F u is an adjustment factor for the forest effect on 
windspeed, e is the vapor pressure, c p is the specific heat, L s is the latent heat of 
sublimation, p a is the air density, and 5y is a stability adjustment. 


Appendix B 

As seen in Section 4B, the mathematical description of the DMRT involves integro- 
differential equations that must be solved numerically and an inverse relationship that is 
difficult to estimate. There are several methods available to determine the inverse. We 
have chosen to use a neural network representation of the DMRT to solve the inverse 
problem. The input-output pairs of the DMRT are used to train the neural network . Once 
the neural network is trained, the brightness temperatures can be computed readily from the 
input parameters S. This has been used for the DMRT in Tsang et al. (1992) and Davis et 
al. (1993). We use a multilayer perceptron (MLP) which is a feedforward neural network 
having one or more layers of hidden neurons between the input and output layers. They 
have a simple layer structure, where successive layers of neurons are fully interconnected, 
with the connecting weights controlling the strength of the connections. Based on the 
snowpack parameters S as the inputs and brightness temperatures T b as the outputs, the 
MLP has four separate input and output nodes. We used one hidden layer with five nodes. 
The MLP is a representation of the function <p where T b = 0(S,F) as shown in Figure 12. 
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An MLP is a function determined by its architecture and its connection weight 
values. The input vector passes directly to the first layer of connection weights. The input 
«,(/+ 1) to each neuron in the next layer is the sum of all of its incoming connection 
weights multiplied by their connecting input neural activation value aj{l). The trainable 
offset value associated with the neuron is added to the sum, and the result is fed into the 
nonlinear function of the neuron, which is usually the sigmoid function (Lippmann, 1987), 

/( W ) = (l + e- u )'\ (B.l) 

The training phase of the MLP uses the backpropagation (BP) learning rule, an 
iterative gradient descent algorithm designed to minimize the mean squared error between 
the desired targets and the actual output vectors. The weights are updated until a minimum 
error criterion is reached between the target values T b and the predicted outputs from the 
MLP. 

To provide a training set for the neural network, the DMRT is run for many 
combinations of the four snowpack parameters, discretizing the space appropriately. The 
neural network is trained on this model output set until the overall improvement at each new 
iteration becomes small. 

This provides the forward model which outputs the 1 km pixel brightness 
temperatures, given the snow hydrology parameters. The adjusted brightness 
temperatures, based on the simulated and measured values are used as inputs to the inverted 
neural network, which then gives adjusted snowpack parameters for each 1 km pixel. The 
snow hydrology model then continues onward with these updated parameters, being 
updated for each day in which there are SSM/I observations available. A least mean square 
error scheme is used to update the snow parameters for each pixel using the neural network 
inversion of the DMRT. The snow hydrology model then continues forward with the 
updated parameters until the next SSM/I observation. This is the forward iterative 
inversion neural network technique discussed in Davis et al. (1993). 
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Table 1. Data inputs used and from where they were obtained. 


Data Type 

Agency 

Address 

a) Hydrology model inputs 



Digital Elevation Model file 

U. S. Geological Survey 

Sioux Falls, SD 

Land Use and Land Cover files 

U. S. Geological Survey 

Sioux Falls, SD 

Hourly point meteorological 
measurements 

National Weather Service 

Asheville, NC 

Daily precipitation totals 

National Weather Service 

Asheville, NC 

Basin streamflow measurements 

U. S. Geological Survey 

Sioux Falls, SD 

b) Satellite remote sensing measurements 


SSM/I microwave brightness 
temperatures 

Marshall Space Flight 
Center, NASA 

Huntsville, AL 

c) Ground truth measurements for algorithm validation 


SNOTEL SWE measurements 

Natural Resources 
Conservation Service 

Portland, OR 
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Table 2. Name, location, and elevation of the SNOTEL stations-available for the year of 
study that are in or near the basin (see Figure 2). 


Number 

Name 

Latitude 

Longitude 

Elevation 



(Deg-Min) 

(Deg-Min) 

(m) 

1 

Beartown 

37-43 

107-30 

3530 

2 

Middle Creek 

37-37 

107-02 

3420 

3 

Slumgullion 

37-59 

107-12 

3470 

4 

Upper Rio Grande 

37-43 

107-15 

2860 

5 

Upper San Juan 

37-29 

106-50 

3080 

6 

Wolf Creek Summit 

37-29 

106-48 

3345 


Table 3. The values used for the fixed parameter vector F in the DMRT for SNOTEL 
station #2 and for station #4. 

Number p 

2 

T, 

Q 

Q* 

(cm) 

0.0375 

(3.0 + iOA)e 0 

4 

270. 

0.40 

0.25 

0.0450 

0.0600 

0.0350 

(3.0 + t0.1)£„ 

270. 

0.40 

0.30 

0.0425 

0.0500 
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Figure 1 Flow chart outlining the steps in the combination of the snow hydrology model 
and remote sensing information to estimate snowpack parameters. Letters refer to 
paragraphs in the text. 
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Figure 2 Location of the upper Rio Grande river basin near Del Norte, CO within the 
region. The left border is the boundary between Utah and Colorado, and the lower border 
the boundary between New Mexico and Colorado. 
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Figure 3 Elevation image of the Upper Rio Grande river basin near Del Norte, CO. The 
image resolution is 1 km based on the DEM. The location of the SNOTEL measurement 
stations are also included (see Table 1). 
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Figure 4 Sensor footprint coverage of a SNOTEL site for four consecutive days. The day- 
to-day variation in the coverage position of the footprint for the SNOTEL pixel is shown. 
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Figure 5 Comparison between hydrology model pixel simulated SWE and measured SWE 
for SNOTEL stations #2 and #4. 
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Figure 6a-d Evolution in time at station #2, high elevation station, of the four snowpack 
state variables in S: depth, temperature, grain size, and fractional volume. 
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Figure 7a-d Evolution in time at station #4, low elevation station, of the four snowpack 
state variables in S: depth, temperature, grain size, and fractional volume. 
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Figure 8 Passive remote sensing with observation of brightness temperatures of a layer of 
dielectric particles embedded in a background medium of permittivity e overlaying a 
homogeneous half-space of dielectric medium of permittivity e 2 . 
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Figure 10 Comparison of the SWE simulations at SNOTEL station #4 for the hydrology 
model without updating and the hydrology model with updating. Four brightness 
temperature channels are used in the multi-parametric inversion. The 37 H brightness 
temperature measurements are shown for reference. 
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Figure 1 1 Comparison of the SWE simulations at SNOTEL station #2 for the hydrology 
model without updating and the hydrology model with updating. Four brightness 
temperature channels are used in the multi-parametric inversion. The 37 H brightness 
temperature measurements are shown for reference. 
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Figure 12 The structure of the multilayer perceptron neural network used to represent the 
dense medium radiative transfer model (DMRT). We use four input nodes, four output 
nodes, representing the input and output vectors S and Tj,, and one hidden layer with five 
nodes. 



